load('data/DailyP.RData')
head(daily_p)
## # A tibble: 6 x 8
## date station lon lat mean_temp min_temp max_temp daily_p
## <date> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2018-10-01 04V -106. 38.1 55.7 46.4 71.6 0
## 2 2018-10-01 0CO -106. 39.8 38.5 33.8 46.4 0.04
## 3 2018-10-01 20V -106. 40.1 54.7 39.2 69.8 0
## 4 2018-10-01 2V5 -102. 40.1 50.9 44.6 68 0
## 5 2018-10-01 2V6 -103. 40.1 48.7 42.8 68 0
## 6 2018-10-01 33V -106. 40.8 54.4 37.4 66.2 0
unique_asos <- daily_p %>%
distinct(lon, lat, station) %>%
st_as_sf(., coords = c('lon','lat'), crs = 4326) %>%
get_elev_point()
## Downloading point elevations:
## Note: Elevation units are in meters
st_read('data/unique_asso_elev.gpkg')
## Reading layer `unique_asso_elev' from data source
## `C:\R-4.1.3\ESS580A9\Metals\data\unique_asso_elev.gpkg' using driver `GPKG'
## Simple feature collection with 67 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -108.7593 ymin: 37.1515 xmax: -102.241 ymax: 40.7503
## Geodetic CRS: WGS 84
monthly_p <- daily_p %>%
mutate(month = month(date)) %>%
group_by(month, station) %>%
summarize(monthly_p = sum(daily_p)) %>%
left_join(unique_asos)
## `summarise()` has grouped output by 'month'. You can override using the
## `.groups` argument.
## Joining, by = "station"
ggplot(monthly_p, aes(x = elevation, y = monthly_p, color = month)) +
scale_color_viridis_c() +
geom_point()
monthly_t <- daily_p %>%
mutate(month = month(date)) %>%
group_by(month, station) %>%
dplyr::select(-lon, -lat) %>%
summarize(across(where(is.numeric), mean)) %>%
left_join(unique_asos,.)
## `summarise()` has grouped output by 'month'. You can override using the
## `.groups` argument.
## Joining, by = "station"
ggplot(monthly_t, aes(y = mean_temp, x = elevation, color = month)) +
geom_point() +
scale_color_viridis_c()
## Warning: Removed 31 rows containing missing values (geom_point).
July_p <- daily_p %>%
st_as_sf(., coords = c('lon','lat'), crs = 4326) %>%
mutate(month = month(date)) %>%
filter(month == 7) %>%
group_by(month, station) %>%
summarize(precip = sum(daily_p))
## `summarise()` has grouped output by 'month'. You can override using the
## `.groups` argument.
# Class hint
# get the boundary data of Colorado
CO <- spData::us_states %>%
filter(NAME == 'Colorado')
# make a empty grid
CO_2163 <- st_transform(CO, crs=2163)
CO_2163_box <- st_bbox(CO_2163) %>%
st_as_stars(dx = 1000) %>%
na.omit(.)
# fit the coordinate reference system
July_p <- July_p %>%
st_transform(., st_crs(CO_2163_box)) %>%
na.omit(.)
# interpolation
interp_July_p <- idw(precip ~ 1, July_p, CO_2163_box) %>%
dplyr::select(1)
## [inverse distance weighted interpolation]
# using tmap
tm_shape(interp_July_p) +
tm_raster(palette = 'plasma', style = 'cont',
title="Total precipitation \n(in) in July 2019") +
tm_legend(legend.outside=TRUE)
# using mapview
mapview(July_p, zcol='precip') +
mapview(interp_July_p, na.col=NA,
col.regions = mapviewGetOption("vector.palette"))
Hint! Use get_elev_raster
# get elevation raster and make grid
ras <- get_elev_raster(unique_asos, z = 8) %>%
raster::crop(., unique_asos)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
ras.star <- st_as_stars(ras, dx = 1000)
names(ras.star) <- "elevation"
# join elevation and precipitation data
unique_asos_July <- st_read('data/unique_asso_elev.gpkg')
## Reading layer `unique_asso_elev' from data source
## `C:\R-4.1.3\ESS580A9\Metals\data\unique_asso_elev.gpkg' using driver `GPKG'
## Simple feature collection with 67 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -108.7593 ymin: 37.1515 xmax: -102.241 ymax: 40.7503
## Geodetic CRS: WGS 84
July_pe <- daily_p %>%
mutate(month = month(date)) %>%
group_by(month, station) %>%
filter(month == 7) %>%
summarise(precip = sum(daily_p)) %>%
left_join(unique_asos_July,.) %>%
na.omit(.)
## `summarise()` has grouped output by 'month'. You can override using the
## `.groups` argument.
## Joining, by = "station"
# Interpolation
interp_July_p_elev <- idw(precip ~ elevation, July_pe, ras.star) %>% dplyr::select(1)
## [ordinary or weighted least squares prediction]
## You will need to create a Stars raster that has elevation data.
tm_shape(interp_July_p_elev) +
tm_raster(n = 10, palette = "-plasma", midpoint = TRUE,
title="Total precipitation \n(in) in July 2019") +
tm_legend(legend.outside=TRUE)
## stars object downsampled to 1346 by 743 cells. See tm_shape manual (argument raster.downsample)
How close do our simple approaches come to reproducing prism maps?